NOTICE: this is the authors version of a work that was accepted for publication in Physica A: Statistical Mechanics 
and its Applications. Changes resulting from the publishing process, such as peer review, editing, corrections, 
structural formatting, and other quality control mechanisms may not be reflected in this document. Changes may 
have been made to this work since it was submitted for publication. A definitive version was subsequently published 
in Physica A: Statistical Mechanics and its Applications, [VOL 391, ISSUE 6, (15 March 2012)] 
DOI:10.1016/j.physa.2011. 11.058 

Critical Temperature Studies of the Anisotropic Bi- and Multilayer Heisenberg 

Ferromagnets in Pair Approximation 

Karol Szalowski 3 *, Tadeusz Balcerzak 3 

"Department of Solid State Physics, University of Lodz, III. Pomorska 149/153, 90-236 Lodz, Poland 

en" 

o 

(N 

Mh- The Pair Approximation method is applied to studies of the bilayer and multilayer magnetic systems with simple cubic 

structure. The method allows to take into account quantum effects related with non-Ising couplings. The paper adopts 
r — ■ the anisotropic Heisenberg model for spin 5=1/2 and considers the phase transition temperatures as a function of 

the exchange integrals strength in line with the role of intra- and interplanar anisotropic interactions in the onset of 
low-dimensional magnetism. The compensation effect for the Curie temperature is found for asymmetric interactions 
within the neighbouring planes of the bilayer system. The paper predicts the saturation of the Curie temperature for 
strong interplanar interactions. However, such an effect for the multilayer system occurs only when the interplanar 
interactions are of purely isotropic character. 
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1. Introduction 

Layered magnets, which are composed of some number of two-dimensional planes, have attracted considerable 
attention both theoretical and experimental [I]. These systems bridge the gap between the two- and three-dimensional 
magnets, so that studying their properties provides some insight into the cross-over between the behaviours character- 
istic of various dimensionalities |2]. That makes the systems particularly interesting from the theoretical point of view. 
Furthermore, experimental realizations of the layered magnets may be based, for example, on organic compounds, 
which offer tunable magnetic properties and a huge variety of structures yj, [3|-|5|] . 

^ The up to date available exact solutions for magnetic models are limited to one-dimensional and some two- 

dimensional magnetic systems, which are mainly of the Ising type jg]. The use of some general transformations 
allows to extend the range of soluble models |7[]. However, in order to study the transition between the two- and 

t^J- | three-dimensional systems, the approximate methods are needed. The desired methods should be able to cope with 

the quantum effects resulting from the non-Ising couplings between the spins. Moreover, such methods should be in 
accordance with the Mermin-Wagner theorem J8J, which imposes rather strict limitations on the magnetic ordering of 
low-dimensional magnets. 

Numerous studies have been mainly concentrated on bilayer ferromagnets, which are the first stage between the 
two-dimensional magnetic plane and three-dimensional bulk magnet. Among such systems, the bilayer Ising model 
on the square lattice has attracted considerable interest I19 H16I1 . The studies included also other underlying planar 



lattice s 1 1711 . especially the Bethe lattice 1 1 8142 111 . Some works covered the interesting case of the spins larger that 1/2 



|9l|22j,|23j], as well as the spins different in both magnetic planes M24h26I1 . The considerations have been extended to 



include dilution 1271 I28H . amorphous structure H291I30H . and the case of unequal coupling strengths in each magnetic 
plane 13111 . The layered systems with anisotropic Heisenberg interactions have been studied in Refs.|5|,[32|]. The prop- 
erties of some more complicated multilayer structures were also investigated 133143711 . In addition, layered quantum 
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Figure 1: (a) A schematic view of the bilayer composed of two layers, A and B. The intraplanar couplings are J^ = Jy = J AA , J AA and 
J BB = J BB = J BB , J BB , respectively. The interplanar coupling is J AB = J AB = J AB , J AB . (b) A schematic view of a multilayer, containing an 
infinite number of subsequent layers A and B. 



Heisenberg antiferromagnets attract the attention, especially as they exhibit a phenomenon of quantum criticallity. 
For example, they are studied by means of Quantum Monte Carlo J38J 13911 or renormalization-group methods |]2j. 
However, according to our knowledge, there are only very few works able to describe (and compare) simultaneously 
within one method the single plane, bilayer, multilayer and bulk systems. 

In order to fill the gap our work describes the above systems in a thermodynamically consistent way. However, 
in this paper we have concentrated mainly on the ferromagnetic bilayer and multilayer systems consisting of two 
magnetically distinct kinds of magnetic planes. A dopting the anisotropic quantum Heisenberg model with spin 1/2, 
the description is based on the Pair Approximation I140[ kill method, which has been substantially modified in order to 
account for the anisotropic structure of the bilayer/multilayer system. We focus our study on the critical temperature, 
and emphasize the influence of various coupling strengths and the importance of interaction anisotropies. Let us 
note that the developed approach allows for a complete construction of the termodynamic description based on the 
expression for Gibbs free energy for the system. 

The paper is organized as follows: in the theoretical section (II) we develop the Pair Approximation method for 
the anisotropic layered systems. In the third section (III) we present the results of numerical calculations for the phase 
transition temperature of the bilayer and multilayer systems. A comparison between these results as well as some 
discussion is offered there. In the last section (IV) the final remarks and conclusion are presented. 



2. Theoretical model 

A schematic picture of the bilayer and multilayer system in question is presented in the Fig.l. The bilayer system 
consists of two parallel planes, A and B, containing A and fi-type of atoms, respectively. We assume that the lattice 
sites form the simple cubic structure. Both intra- and interplanar exchange interactions are of the anisotropic Heisen- 
berg type. We consider the case when < f^ = Jl^ - f± < JT (v = A, B; /u — A, B), which covers all intermediate 
situations between the pure Ising {J v ^ - 0) and the isotropic Heisenberg (/* = J* 1 *) models. 

The multilayer system is built from the interacting bilayer segments by the infinite repetition of the bilayer in the 
direction perpendicular to the A and B planes. 

The Hamiltonian of the bilayer can be written in the form of: 



(ieA.jeB) (ieB.jeB) 

- J] {j^isisi + s^s^ + j^sisll-iiYsi-iiYsL 

(ieAJeB) ieA ieB 



(1) 
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where h stands for the external magnetic field. In the case of multilayer system the bilayer Hamiltonian should be 
summed up over all A and B layers, where each layer interacts with its neighbours from both sides. 

The general method of constructing the thermodynamic description within PA has been sketched out in our pre- 
vious papers ikdklll . The crucial point of the PA is the method of construction of the single-site and pair density 
matrices, p' ev and p l€v -J € f (y = A,B; p - A, B), respectively. In the present case the matrices take the form of: 



(v = A, B), and 



piev = ^G" exp y p (A v + h ) S f v ] 



where the pair Hamiltonian is of the form: 

Tffevj*, = jY^sfSi^ + S^Si^ + J^SfSf" 

+ (A v " + A v/J /2 + h) S i& + (A v/J - A v "/2 + h) S 



./&/' 



(2) 
(3) 

(4) 



{v = A,B;n=A,E). 

The single-site Gibbs energies, G A and G B , which are obtained from the normalization condition for the single-site 
density matrices p ,eA andp' eB , respectively, are given by the expressions: 

(5) 

(6) 

On the other hand, the pair Gibbs energies, G AA , G BB and G AB , which are obtained from the normalization condition 
after diagonalization of the pair density matrices p l€A <J €A < p' eB J eB and p l€A J €B t respectively, can be presented as: 



G A = 


-k B T In < 2 cosh 


[!(*+»)] 


G B = 


-k B T In < 2 cosh 


|(A« + »); 



G AA = -k B T In \l exp [^f~\ cosh [/3 (a^ + h)] + 2 exp (-^-) cosh (p^) 
G BB = -k B T In 1 2 exp V-^— J cosh [b (a bb + h)] + 2 exp (-^H cosh (wf 5 ) 



(7) 
(8) 
(9) 



Application of the PA method leads to the expression for the total Gibbs energy, from which all thermodynamic 
properties can self-consistently be derived. The general formula for the Gibbs energy has the form of: 



G AB = -k B T In J 2 exp l^pj cosh [b (a ab + h)] + 2 exp (-^f-) cosh | ^(A AB f + (j AB ) 






Qiev.jen _ Z £ Iq<^v + Qfif\ 



(10) 



(v = A, B; p — A, B), where G' ev = G y are the single-site and G' ev/e/J = G v,/J are the pair Gibbs energies, respectively. 
The number of nearest neighbours (NN) is denoted by z (z - 5 for the bilayer and z — 6 for the multilayer structure). 
Taking into account the total number of NN pairs, the Gibbs energy per 1 lattice site for the bilayer system is given 
by: 

G = G AA +G BB + G AB /2-2(G A +G B ) (11) 

whereas for the multilayer system we obtain: 



G = G M + G BB + G AB - -(G A + G B ). 



(12) 



The density matrices, and thus the Gibbs energy, depend on the parameters A v , A VjU and A VyU , which can be 
expressed by means of four other parameters X ,Ji . The determination of these parameters is presented in full detail in 
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the Appendix A. In principle, they are calculated from the variational minimization of the total Gibbs energy, which 
yields the conditions Tr (p' ev S' : ev ) - Tr [p' ev -^S l ey \ Having obtained the parameters, the Gibbs energies ([TO and 
(fT2b and hence the total Gibbs energy ( ITOb are at our disposal as the functions of arbitrary temperature (B - l/k B T) 
and the external magnetic field h. Let us emphasize that the local anisotropic exchange interactions (f^ and J v f) 
are taken into account in the Eqs.( lA.4t - (!A.7b determining the parameters A v '^. Thus, these local equations give a 
possibility to study a variety of interesting cases for the intra- and interplanar couplings. 

The main focus of the present work is the critical temperature T c of the bilayer and multilayer systems. The 
procedure used to determine T c (in case of the second-order phase transition) for the systems of interest is presented 
in the Appendix B. In a general case, calculation of T c requires solving a determinant equation (IB.U of the size 4x4. 
However, in the particular case when the two layers A and B are magnetically equivalent, i.e., for J^ - J BB - J z , 
J AA - J B B - J ± ,we obtain the simple formulas for the critical temperature: 

4C + C 4B = 3 (13) 

for the bilayer system, and 

2C + C 4B = 2 (14) 

for the multilayer system. The C and C AB coefficients take the form of: 

J± 



C - exp -— — — cosh., 

y \ 2k B T c J \2k B T c 

I jAB \ I jAB \ 

On the basis of Eqs.(fT3l and ( TBI it is a simple task to check that, for instance, for the uncoupled Ising layers (when 
J± - and J AB - J AB - 0) we obtain: k B T c /J z = 1/2 In 2. On the other hand, when the Ising planes are coupled, 
the result is k B T c /J z = 1/2 ln[z/(z - 2)] where z - 5 is the NN number for the bilayer system, and z — 6 stands for 
the multilayer system. We see that the result for the uncoupled planes corresponds to the Curie temperature of one 
monolayer with s. q. st ructure (z - 4). These analytical solutions are in agreement with the findings of previous papers 
on the PA method H El . 



In other limiting case, when the system contains only isotropic Heisenberg interactions (J ± = J z = J), and 
assembles uncoupled planes, we obtain k B T c /J z = which is in agreement with the Mermin-Wagner theorem J8J. On 
the other hand, for the interacting Heisenberg planes we obtain k B T c IJ : — 1 / ln[z/(z - 4)], where z — 5 stands for the 
bilayer and z — 6 corresponds to the multilayer system. Again, this result is in agreement with previous studies of the 
Heisenberg systems within the PA method 14211 . 

From the above remarks one can conclude that the PA method should be able to reveal the onset of magnetism 
in the case of the isotropic Heisenberg planes (J^ = J BB = J AA = J BB = J) when the interplanar interactions are 
gradually being switched on. It turns out that either from Eq.(fT3l) or (TPfl i. when J AB — > 0, we find the common result: 

(16) 



ln(j? B /j) 



It is interesting to see that in the logarithmic law (fTol l the perpendicular interaction between the planes {J AB ) has no 

influence on the Curie temperature. 

Finally, for the non-interacting planes (J AB = J AB - 0), when we define the intraplanar anisotropy parameter as: 

6 - (J, - J±)/J~, we find that: 

k B T c 1 

-5-^ = (17) 

J z In (1/6) V ; 

for 6 — > 0. This result is in agreement with our previous findings for the anisotropic Heisenberg model in the diluted 



systems 14 1H . when the effective coordination number amounts Zes = 4 . It describes the onset of 2D magnetism as a 
result of anisotropic interactions. 

The numerical results based on the general equation for the critical temperature (IB. IK especially where the ana- 
lytical solutions do not exist, are presented in the next section. 
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Figure 2: The critical temperature of a bilayer system with Ising interlayer and intralayer couplings as a function of relative strength of interlayer 
interaction. Dashed line is the MFA result, solid line - PA result. Triangles denote the results of Monte Carlo method (Ref. 11211 ). circles - Comer 
Transfer Matrix Renormalization Group (Ref. HllO . squares - Transfer Matrix Mean Field Approximation (Ref. I13l0 . 



5- 



4 
< 

5 2 

J? 



1 - 



0- 



T ' ■ ■ 1 ■ 

Bilayer Multilayer 



JAB/JAB = : 

J AB /J AB = 




j m /j z m = J BB /J? B = , Jf»/J/* = 1 



8 



12 



16 



Figure 3: The critical temperature of the bilayer and multilayer system with equal Ising couplings within the magnetic planes as a function of the 
relative strength of interlayer interaction. The interlayer interaction is either of the Ising or isotropic Heisenberg type. 



3. Numerical results and discussion 



First let us study the case of the critical temperature of the bilayer system composed of two magnetic planes with 
interplanar couplings of the Ising type and equal exchange integral strength (J BB /J AA = 1, J BB /J? B = J^/J^ = 0). 
The critical temperature is plotted in the Fig. [2] vs. J AB /J AA . The linear dependence predicted by the Mean Field 
Approximation, k^T^ FA - J AA + J AB /4, is shown by the dashed line, while the solid line presents the result of 



the present (PA) method. The closed symbols depict the predictions of Monte Carlo simulations (after Ref. 111211 '). 



le pre 
O) 



Transfer Matrix Mean Field Approximation (after Ref. II 1 311 > and Corner Transfer Matrix Renormalization Group 
(after Ref. Ill ill ). It is visible that the MFA prediction significantly overestimates T c , while the PA results, which form 
a nonlinear dependence, are much closer to the above mentioned accurate estimations of the critical temperature. 

An influence of the increasing interplanar coupling on the critical temperature of the bilayer or multilayer system 
is illustrated in the Fig. [3] The plots show the dependence of k^Tc/J^ on J AB /J AA . It is assumed that the intraplanar 
interactions within A and B planes are of purely Ising type {J AA = J BB - 0) and they are taken with equal strength 
(J AA = J BB ). The AB couplings are assumed to be either of Ising or isotropic Heisenberg type. All the curves 
start from a common point k^T c /J AA - 1/2 In 2, which constitutes the critical temperature of a single Ising plane. 
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Figure 4: The critical temperature of the bilayer and multilayer system with equal and unequal Ising couplings within the magnetic planes as a 
function of the relative strength of interlayer interaction. The interlayer interaction is of the isotropic Heisenberg type. 




Figure 5: The critical temperature for a bilayer system with unequal strength of intralayer couplings in both planes, either for the Ising or isotropic 
Heisenberg coupling between the planes. 
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Both for bi- and multilayer, the Ising coupling elevates the Curie temperature stronger that the isotropic Heisenberg 
coupling. For the case of a bilayer system, the critical temperature tends to saturate even though the coupling J AB jJ AA 
is increased to infinity. The saturation value for the Ising coupling is found as k^TdJ^ = l/21n(4/3), whereas the 
corresponding limiting value for strong isotropic Heisenberg coupling is ksT c /J AA = 1/2 In (8/5). The tendency for 
the critical temperature to saturate is also clear for the case of the multilayer system with the isotropic Heisenberg 
AB interactions. Interestingly, the limiting critical temperature of a strongly Ising-coupled bilayer system is the same 
as for the strongly Heisenberg-coupled multilayer system. On the contrary, for the multilayer system with pure Ising 
coupling the Curie temperature tends to increase unlimitedly when J AB j] AA is increased. Such behaviour of the curves 
is in agreement with Eqs.(fl3T> and (fT4l >. 

The same kind of behaviour has been detected for the case of a bilayer and multilayer system for which the 
intraplanar interactions are of isotropic Heisenberg type . For this reason such temperatures have not been illustrated in 
a separate plot. The only qualitative difference is that all the curves commence at T c equal to zero for uncoupled planes. 
The unlimited increase of the critical temperature takes place only for the multilayer system with Ising coupling. 
Also, the critical temperatures of the bilayer system with the Ising interlayer coupling and for the Heisenberg-coupled 
multilayer system tend to a common limit, as for the previously described case. 

The influence of unequal strength of the Ising intraplanar coupling within the layers A and B on the behaviour of 
bilayer and multilayer system is presented in the Fig.|4] In such asymmetric situation the Curie temperatures have to be 
obtained from the full determinant (IB. It . It is visible that in all the cases considered in Fig.|4]the critical temperature 
tends to saturate for strong interlayer coupling. However, in the case of the bilayer with J BB /J AA =0.1 and strong 
coupling the limiting critical temperature is significantly lower than T c for the single A plane with Ising interactions. 
This suggests that in such a case the presence of the Heisenberg lateral bounds weakens the ferromagnetic order in 
the planar system. An analogous effect is absent in the multilayer system with J BB /J AA =0.1. Moreover, for equal 
strength of intraplanar couplings (J BB /J AA = 1), the critical temperature of multilayer system can be elevated as a 
result of the presence of Heisenberg isotropic coupling between the magnetic planes. 

To emphasize the influence of the isotropic Heisenberg coupling on the critical temperature for a bilayer system 
with unequal intraplanar couplings, we present the next Fig. [5] for J BB /J AA =0.1. It can be observed that in the 
case of the Ising coupling within the planes (J^/J^ - J BB /J BB = 0), the Ising interplanar coupling causes the 
critical temperature to increase, while the isotropic Heisenberg coupling reduces T c . The critical temperature remains 
reduced even for a strong coupling. The situation is qualitatively different for the case of the bilayer system with 
isotropic Heisenberg interactions within the planes (J AA /J AA = J BB /J BB = 0). For the uncoupled system, the critical 
temperature is equal to zero. Initially, introducing either the Ising or isotropic Heisenberg coupling T c becomes 
elevated. For the Ising AB interaction, the temperature keeps increasing monotonically and tends to saturate for 
strong coupling. On the contrary, for the Heisenberg interplanar coupling, T c reaches a certain maximum value and 
then becomes reduced until it reaches some limiting value (quite low when compared to the case of the Ising coupling). 

The behaviour of the critical temperature of the bilayer and multilayer system with isotropic Heisenberg intrapla- 
nar couplings deserves particular attention. As stated in the theoretical part a single plane with isotropic Heisenberg 



interactions does not show magnetic ordering according to the Mermin-Wagner theorem []8j |43J. Therefore, the rise 
of ordered state as a consequence of introducing the interaction between the planes appears especially interesting. In 
the Fig. [6] the critical temperature k^T c /J AA is shown as a function of the inverse logarithm of the couplings ratio 
J AB /J AA . Both the case of bilayer and multilayer system is taken into account. Moreover, the AB coupling is assumed 
either in the Ising or isotropic Heisenberg form. It is visible that for both systems and both kinds of interplanar in- 
teractions the critical temperature follows the same asymptotic behaviour when J AB /J AA vanishes. It has been shown 
analytically in the previous section that this asymptotic behaviour takes the form of k-aT c /J AA = -1/ In (j AB /J AA ). 
For stronger coupling it is visible that the interplanar interaction has a more pronounced effect on the multilayer than 
on the bilayer system. Moreover, the Ising coupling causes the critical temperature to increase more rapidly than the 
isotropic Heisenberg one. Let us note that such kind of behaviour predicted here by the PA method is consistent with 
the expectations based on the scaling theory 114411 as well as with some other approaches ]2j, |45[ |46|] . It bears some 



resemblance to the result for a single Heisenberg plane, in which the interaction anisotropy in spin space breaks the 
assumptions of the Mermin-Wagner theorem, and for which case the PA method predicts T c <x -1/ In [(J z - J ± ) /J z ] 
(see Eq.dTTb and Ref. J41fl ). This result is in agreement with Ref. yfl. 

The asymptotic behaviour of the critical temperature can also be tested for the case of unequal couplings in both 



NOTICE: this is the authors version of a work that was accepted for publication in Physica A: Statistical Mechanics 
and its Applications. Changes resulting from the publishing process, such as peer review, editing, corrections, 
structural formatting, and other quality control mechanisms may not be reflected in this document. Changes may 
have been made to this work since it was submitted for publication. A definitive version was subsequently published 
in Physica A: Statistical Mechanics and its Applications, [VOL 391, ISSUE 6, (15 March 2012)] 
DOLT0.1016/j.physa.2011. 11.058 



0.8- 



i — ' — ' — ' — i — ■ — ' — ' — r 
. Bilayer Multilayer 




0.4 0.6 

-1/ln(J AB /J AA ) 



1.0 



Figure 6: The critical temperature for the bilayer and multilayer system with isotropic Heisenberg intralayer coupling as a function of inverse 
logarithm of the relative interplanar exchange interaction. The interplanar coupling is either of Ising or isotropic Heisenberg type. 
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Figure 7: The critical temperature for a bilayer system with isotropic Heisenberg intralayer coupling, as a function of inverse logarithm of the 
relative interplanar coupling energy. The interplanar coupling is of isotropic Heisenberg type. Either equal or unequal intralayer couplings are 
assumed. 
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Figure 8: The normalized critical temperature for a bilayer system as a function of the relative strength of the Ising couplings within A and B planes 
for various interplanar coupling energies: (a) Ising coupling between the planes; (b) isotropic Heisenberg coupling between the planes. 



layers A and B with isotropic Heisenberg-type interactions. Such a situation for the bilayer system is presented in the 
Fig-ll] It follows that the proportionality of T c to 1/ \n(j AB /J AA ) is unchanged by weakening the relative strength of 
coupling within the B layer; however, the coefficient of this proportionality decreases. 

The case of unequal strength of the intraplanar AA and BB couplings may have also other important consequences. 
For instance, the behaviour of the critical temperature as a function of the ratio J BB /J AA for various interplanar AB 
couplings is especially worth mentioning. Let us fix the value of J AA and normalize the critical temperature to it. In 
the Fig. [Hfa) the dependence of the critical temperature on the ratio of intraplanar couplings J BB /J AA is presented for 
selected values of the pure Ising-like interaction between the planes, i.e. when J AB = 0. It is visible that within the 
limit of vanishing intraplanar couplings J BB the critical temperature of the system tends to the value predicted for 
the sole A magnetic plane, i.e. k^TdJ^ = 1/2 In 2, independently on the energy of the interplanar interaction. The 
increase of J BB elevates T c and the effect is more rapid for stronger coupling between the planes. 

However, the situation noticeably changes when the pure Ising interplanar coupling is replaced by the isotropic 
Heisenberg one (i.e., for J AB = J AB ), which is depicted in the Fig.|8](b). For a very weak or even vanishing coupling in 
the plane B, the effect of Heisenberg coupling J AB on the critical temperature is detrimental; a stronger J AB coupling 
reduces significantly the critical temperature below the value of k^ T C /J AA - 1 /2 In 2 which is characteristic of a single 
A plane. Such a decrease can be interpreted as an effect of the coupling between the perpendicular components of 
the spins, tending to destroy the spin order along the z axis, and is indeed a quantum effect. This demonstrates that, 
in certain circumstances, the ferromagnetic coupling is able to reduce the value of the critical temperature. As it is 
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Figure 9: The position of the compensation point for the bilayer system with the Ising coupling within the planes as a function of the relative 
strength of isotropic Heisenberg interlayer interaction. 



visible in Fig.[8](b) this reduction of T c becomes less pronounced for the increasing coupling in the B plane. On the 
other hand, when J BB is comparable with J AA the behaviour of the critical temperature follows the trend presented in 
the Fig.[8la). It is particularly interesting that at some value of J BB jJ AA the critical temperature k^T c /J AA - 1/2 ln2 
is restored for each curve, even though the disordering Heisenberg interplanar couplings are present. The value of 
J BB /J AA at which the compensation of the opposite influences of J AB and J BB takes place is very weakly sensitive to 
the strength of interplanar coupling. This is visible in the Fig. |8](b) as an approximate interception point of all the 
curves plotted for various values of J AB /J AA . 

The insensitivity of the coupling ratio (j BB /J AA ) , for which T c has the constant value (characteristic of an uncou- 
pled bilayer system), on the interplanar interaction is further illustrated in the Fig. [9] There, the values of (j BB /J AA ) 
are plotted as a function of the interplanar coupling J AB /J AA strength. The plot is prepared not only for isotropic 
Heisenberg interplanar coupling (J AB jJ AB — 1), but also for two smaller values of interplanar interaction J AB . It is 
noticeable that the curves are almost flat (which is most evident for the isotropic Heisenberg AB interaction) so that 
the compensation point U BB /J AA ) is quite well defined for a wide range of J AB . The plot also shows that the value of 
U BB /J AA ) is shifted towards lower values when J AB /J AB decreases. In the limiting case of Ising AB coupling, i.e., 
for J AB = 0, the 'compensation' takes place at J BB /J AA - 0, which is just the result illustrated in the Fig. 0a). 

Finally, let us study how the relative amount of perpendicular component J AB in the AB couplings influences the 
positional diffusion of the compensation point. Fig. [l0]presents the values of U BB /J AA ) for J AB /J AB varying from 1 
(isotropic Heisenberg AB coupling) down to (Ising coupling). For each particular value of J AB /J AB all the results 

U BB /J AA ) , which are adequate for J AB /J AA ranging from to 1 (i.e. covering the range of the Fig. |9}, are shown by 
the shadowed area. The 'thickness' of this shadowed area plotted in the Fig.[10]corresponds then to the variation of 
(j BB /J AA ) as a function of J AB IJ AA . Thus, the smallest thickness means the weakest sensitivity of the compensation 
point to the strength of interplanar coupling. It is visible that the compensation points are only slightly sensitive to AB 
coupling strength either for the isotropic Heisenberg interaction or when J AB — > 0, while for the intermediate values 
of J AB /J AB the compensation is considerably diffused. 

In the same picture (Fig. ITOb we also included the case when the intraplanar interactions within A and B planes 
are not purely of the Ising type. It is observable that a similar compensation takes also place when J AA /J AA - 

J BB /J BB > 0, but the effect is less pronounced and the value of U BB /J AA ) needed to compensate the influence of AB 
coupling on the critical temperature is lower. We have observed that the critical temperature of a bilayer ferromagnetic 
system undergoes the strongest reduction when the interplanar coupling is of the isotropic Heisenberg type and the 
intraplanar coupling within the planes is of the Ising type. Let us also mention that the compensation phenomenon is 
characteristic of a bilayer system with J AB /J AB > 0, and no similar effect has been observed for a multilayer system, 
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Figure 10: The diffusion of the compensation point for a bilayer system as a function of the interlayer interaction anisotropy. The intraplanar 
couplings were assumed either as the Ising or anisotropic Heisenberg ones. 



when the system consists of an infinite number of subsequent A and B planes. 



4. Final remarks and conclusion 

It has been found that the PA method is useful for the magnetic studies of both the bilayer and multilayer sys- 
tems. In particular, the ferromagnetic phase transition temperature has been thoroughly examined. The anisotropic 
Heisenberg model has been adopted with the variable exchange integral between the planes. 

The model allows to control two limiting cases of the system: the first case is when the planes are totally separated 
and their properties are purely 2-dimensional. The second limiting case comprises situation when the system becomes 
an infinite bulk crystal. The application of the PA method for those limits gives a good consensus with the results 
existing in the literature. In particular, for 2D Heisenberg systems the agreement with Mermin- Wagner theorem has 
been found. 

For an intermediate situation when the bilayer or multilayer system is considered, the PA method allows to take 
into account not only the magnetic anisotropy in spin space (with various J ± and /.), but also the directional anisotropy 
(different J AA , J BB and J AB ) as well. In particular, it enables to study the layers with asymmetric exchange interactions 
(J AA ± J BB ). Thus, the current approach is formulated more generally than that presented in the previous work 
(Ref. JU). 

Regarding the main results obtained in the paper, the asymptotic law for the Curie temperature decay (Eq. Q~6), 
when the interplanar coupling decreases, has been derived. It has been found that only the Ising-type coupling influ- 
ences the Curie temperature in such asymptotic law, and the formula is common both for the bilayer and multilayer 
system. This result of the PA method is in agreement with the predictions of the scaling theory IH411 as well as with 
some other approaches 1121 |45hI47 I . 

Another important result concerns the asymptotic behaviour of the Curie temperature when the interplanar cou- 
pling increases to infinity. It has been found that for the bilayer with intraplanar interactions both of the Ising and 
Heisenberg-type, when the interplanar (either Ising or Heisenberg) coupling increases, the Curie temperature reaches 
its upper limit. The limiting values can be found on the basis of Eq.(fT3b. On the other hand, for the multilayer sys- 
tem with intraplanar interactions both of the Ising and Heisenberg-type, the saturation of the Curie temperature takes 
place only when the interplanar couplings are of the Heisenberg-type. In this case, the limiting values of the Curie 
temperature can be found from Eq.(fT4t. 



Finally, for the bilayer system with asymmetric intraplanar interactions (/ 



J ) the 'compensation point' 



for the Curie temperature is found, as shown in the Figj8](b). This effect arises from some decreasing of the Curie 
temperature when the / BB -interaction decreases, provided the interplanar couplings are of the Heisenberg type. Such 
an effect has not been found in the multilayer system. It reflects the fact that the properties of multilayer system cannot 
be considered as a simple multiplication of the bilayer properties, but they can be substantially different. In the view 
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of that finding a confirmation of the 'compensation point' existence by other methods, for instance, by the quantum 
MC simulations would be desirable. 

Let us remark that according to the theoretical part the application of the PA method far exceeds the calculations of 
the phase transition point. Since the Gibbs energy is at our disposal, the calculations of all thermodynamic quantities 
for the system in question is possible. For instance, the studies of the thermodynamic response functions, such as the 
susceptibility or specific heat might be of interest. However, in our opinion, such studies can be done in a separate 
paper. 

The method can be developed to include antiferromagnetic interactions. Further extensions may also comprise the 
structurally disordered systems as well as the models containing higher spins. 
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Appendix A. Determination of variational parameters 

The molecular field parameters A and A are specified separately for the bilayer and multilayer system, taking into 
account the different number of NN bounds. For instance, for the bilayer system: 

A A = 4A AA +A AB 

A B = 4A BB + A BA 

A AA = 3A AA +A AB 

A BB = 3A BB + A BA 

A AB = 2 (X AA +A BB ) 

A AB = 4 ( A AA_ A BB'j 

A AA = A BB=q (A1) 

In turn, for the multilayer system the parameters are of the form: 

A A = 4A AA +2A AB 

A B = 4A BB + 2A BA 

A AA _ 3A AA + 2A AB 
A BB = 3A BB + 2A BA 

A AB = 2 ( A AA +A BB} + l(jAB +A BAj 

A AB = A ( A AA_ A BB^ + A AB_ A BA 

A AA = A BB = Q (A2) 

The four parameters A w (v = A, B; yu = A, B) introduced in Eqs. dA.U and dA.2b have a meaning of the one-bond 
molecular fields, i.e., the fields acting on the v-atom and originating from its yu-neighbour. These parameters can be 
obtained from the variational principle for the free energy : 

dG 

= (A.3) 



dA'v 
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Let us remark that in the PA method the variational principle [A3]is equivalent to the condition that the mean value 
of a given spin (S l ev ) (magnetization) remains the same when it is calculated either from p' ev or p' €v 'J € f density matrix 

i.e., Tr (p ,ev S! ev J = Tr (p' ev/e/J 5'i ev J. We emphasize that this condition is vital to the correct description of the system, 
since it guarantees the consistence between the quantum state of a single lattice site described either by single-site or 
by pair density matrix. Let us mention that for spin 1/2 such condition guarantees also the reduction property of the 
density matrices p' £v - Tr, (p iev '^\ 149J. Thus, from Eq.( lA.3l l we obtain the set of four variational equations: 



tanh 



tanh 



tanh 



tanh 



f(A A + *) 



?K+*) 



SK + *) 



? (*'+*) 



ex P (^r) sinh [j6 { aAA + h)] 



exp(^)cosh \p(A AA + It)] + exp(-^)cosh(^) 



exp(^)sinh[/?(A AB + /j)] 



V(A«) 2 + (yf) 



= exp 



(-'-?) 



sinh 



W( A ") 2+ w) 2, l 



cosh 



exp (&£■) cosh [p (A AB + It)] + exp (-^f\ 

sxp(^f-)smh[/3(A BB + h)] 
exp (&ZL \ cosh \j3 (A BB + h)] + exp(-^)cosh(^) 



J 8 1 /(A' iB ) 2 +(-'f) : 



exp (^) sinh [/?(A AB + /i)] 



V(A-"0 2 + (/f) 2 



exp 



(-'-?) 



sinh 



i^f+WY 



exp (&tf-) cosh [J3 (A AB + h)] + exp (-^f-) cosh 



/iV(A''«) 2 + (yf) 2 ' 



(A.4) 



(A.5) 



(A.6) 



(A.7) 



The set of Eqs. (lA.4t - (!A.7t has the same form both for the bilayer and multilayer system; however, the difference 
lies in the molecular field parameters A and A which are different for those two cases (and are given by Eq. dA.lt 
or Eq. (IA.2b . respectively). From the above variational equations the parameters A AA , A BB , A AB and A BA can be found 
numerically. 



Appendix B. Determination of the critical temperature 

One of the straightforward applications of variational equations (!A.4t - (!A.7b is the critical temperature calculation. 
For the continuous phase transitions the molecular field parameters should vanish. As a result of linearization of 
Eqs. (lA.4b - (!A.7b (when A w — > 0) the critical temperature T c can be determined from the condition: 

det(M) = (B.l) 

where M is the characteristic matrix (4 x 4) whose elements for the bilayer system can be written as: 

exp [p c J AA 12) - 2 cosh (j3 c J AA /2) 



M, 



Mv 



= M 34 =fl 



= M 33 = 



exp (/3 C J AA I2) + cosh (^7^/2) 
Pc exp(p c J AA /2) - cosh(/Wf 12) 



M 



13 - 



M 2 i - M 44 - 



2 exp (p c J AA /2) + coshfaj^ /2) 
M14 = M23 = Mn - M i2 - M42 - 

(4/J AB ) sinh (p c J AB /2) - 2p c cosh (p c J AB : /2) 



exp (p c J AB /2) + cosh (PcJi 8 / 2 ) 



13 



NOTICE: this is the authors version of a work that was accepted for publication in Physica A: Statistical Mechanics 
and its Applications. Changes resulting from the publishing process, such as peer review, editing, corrections, 
structural formatting, and other quality control mechanisms may not be reflected in this document. Changes may 
have been made to this work since it was submitted for publication. A definitive version was subsequently published 
in Physica A: Statistical Mechanics and its Applications, [VOL 391, ISSUE 6, (15 March 2012)] 
DOIT0.1016/j.physa.2011. 11.058 

M 22 = M 43 = -y 

_ 2B C exp (b c J a ± b /2) - (4/Ji B ) sinh (p c J^/2) 

exp (b c J ab /2) + cosh (/3 c J AB /2) 

(B.2) 



where B c = l/k B T c . 

For the multilayer system the matrix elements take the form of: 

exp (j3 c J AA /2) - 2 cosh (b c J aa /2) 



M l2 = B, 



exp(p c J AA /2) + cosh (pj^ll) 
exp (b c J aa /2) - cosh (/3 c J AA /2) 



exp (j3 c J AA /2) + cosh (/3 c /^/2) 
Mn = M14 - Mn - Mi2 - 

(4/J AB ) sinh (/3 c J AB /2) - 2B C cosh (b c J ab /2) 



M 2 i = M 44 

M 2 2 = M43 

M 23 = M 42 



exp (b c J ab /2) + cosh (/3 c Ji B /2) 
(l/J AB )smh(/3 c Ji B /2) /3 c cosh(/3 c J AB /2) + (J3 c /2)exp(j3 c jf B /2) 



M 24 = M41 = 



exp {p c J AB l2) + cosh (j3 c J AB /2) exp (f3 c J AB /l) + cosh (j3 c J AB /2) 

(J3 c /2)exp({3 c J AB /2) - (l/J AB )smh(/3 c J AB /2) 

exp (f3 c J AB /2) + cosh {/3 C J AB I2) 
2J3 C exp (/3 c J AB /2) - (4//^) sinh (j3 c J AB /2) 



M 33 = A 



exp (f3cJ AB /2) + cosh (/3 c J AB /2) 
exp (/? e 7f B /2) - cosh (/Wf /2) 
exp (f3 c J BB /2) + cosh (f3cJ BB /2) 



exp (p c J BB /2) - 2 cosh (/3 c J BB /2) 

M 34 = & V- ^r / RR / (B.3) 

exp[/3 c J? B /2) + cosh (/3 c J BB /2) 

In a general case the determinant equation (1B.U can be solved only numerically. In the particular case, when the 
layers A and B are magnetically equivalent, the determinant can be reduced to the form of eqs. ([T"3T i or (TT4l . It should 
also be noted that the physical solution corresponding to the critical temperature of the system is the highest value of 
T c obtained from ( IB. lb . 
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